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Contrary to standard coupled cluster doubles (CCD) and Brueckner doubles (BD), singlet-paired 
analogues of CCD and BD (denoted here as CCDO and BDO) do not break down when static corre¬ 
lation is present, but neglect substantial amounts of dynamic correlation. In fact, CCDO and BDO 
do not account for any contributions from multielectron excitations involving only same-spin elec¬ 
trons at all. We exploit this feature to add—without introducing double counting, self-interaction, 
or increase in cost—the missing correlation to these methods via meta-GGA density functionals 
(TPSS and SCAN). Furthermore, we improve upon these CCDO-I-DFT blends by invoking range 
separation: the short- and long-range correlations absent in CCDO/BDO are evaluated with DFT 
and the direct random phase approximation (dRPA), respectively. This corrects the description 
of long-range van der Waals forces. Comprehensive benchmarking shows that the combinations 
presented here are very accurate for weakly correlated systems, while also providing a reasonable 
description of strongly correlated problems without resorting to symmetry breaking. 


I. INTRODUCTION 

Although coupled cluster (CC) theory is one of the pil¬ 
lars of quantum chemistry, commonly used CC methods 
{e.g., singles and doubles, CCSD) often fail in the pres¬ 
ence of static (or strong) correlation [I] [5] . These failures 
are connected to instabilities that appear in the random 
phase approximation (RPA) when Hartree-Fock (HE) 
becomes unstable towards symmetry breaking [5]— RPA- 
like terms are contained in the equations of traditional 
coupled cluster techniques [SHS] • (RPA instabilities come 
in two flavors: ph-ph and pp-hh. The former is related to 
spin and the latter to number symmetry breaking. “In¬ 
stability” in the CC context means that the solution to 
the CC equations results in large or complex cluster am¬ 
plitudes and unbound correlation energies.) Such insta¬ 
bilities may be avoided by modifying the cluster operator 
to include only singlet-paired operators as in, e.g., pair 
CCD [6UTT] (pCCD) and the recently proposed CCDO 
method [5]. However, this workaround is not a panacea: 
the suppression of terms in the cluster operator inevitably 
leads to neglecting part of the correlation. For example, 
CCDO contains none of the correlation terms arising from 
excitations involving same-spin electron pairs which are 
present in CCD (see below). Despite this clear problem, 
pCCD and CCDO are proper wavefunction theories: both 
the ansatz and its solution are rigorous approximations 
to solving Schrodinger equation. The main objective of 
the present work is to investigate ways to add the corre¬ 
lation that is missing in CCDO via physically-motivated 
combinations of this method, and its extension using 
Brueckner orbitals (denoted as BDO), with techniques 
based on density functional theory (DFT) and the RPA. 

The idea of merging techniques like CCDO with DFT 


is inspired by a large number of previous works |121445] 
in which multireference (MR) methods {e.g., complete 
active space or CAS) are mixed with density functional 
approximations (DFAs), in order to account for static 
and dynamic correlation simultaneously. The motiva¬ 
tion for MR-I-DFT is that the MR ansdtze typically used 
in quantum chemistry capture mostly (albeit not exclu¬ 
sively) static correlation, whereas common DFAs provide 
an efficient alternative to evaluate the dynamic correla¬ 
tion. Under this premise, one could think of approximat¬ 
ing the exact energy as 

E + Ef^^[n], (1) 

where E^^ is the MR energy and E^^^ is a (dynamic) 
correlation functional of the density n. Similarly, singlet- 
paired coupled cluster methods (CCO) can describe static 
correlation [5] and hence the incentive for CCO-t-DFT is 
analogous to that of MR-l-DFT. Nevertheless, despite the 
apparent simplicity and soundness of MR-fDFT, these 
methods have not achieved widespread use even though 
concrete implementations have appeared in the literature 
since circa 1970 [121 [13]. This is largely because of three 
problems that plague MR-I-DFT: 

1. Double counting. This problem arises because, in 
general, the MR energy contains some dynamical 
correlation which is also described by the DFA. A 
raw implementation of Eq. [fusing a standard DFA 
will therefore often yield too low energies and un¬ 
satisfactory results [25l [43]. Strategies that have 
been proposed to avoid the double counting include 
the use of small active spaces and reparametrized 
functionals [T^14j : scaling the correlation energy 
density by factors based on the local density ap¬ 
proximation P51 - I5D] (LDA) or on the cusp con- 
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dition for the pair density Ml partition of the 
electron-electron interaction to create MR-I-DFT 
global or range-separated hybrids; 

separation of the correlation contributions to the 
single-particle spectrum inilSl]; as well partitions 
of the orbital space [44l |45] based on embedding 
schemes [M1I17]- None of these solutions is per¬ 
fect: they either eliminate double counting approx¬ 
imately only, or eliminate it exactly but may ne¬ 
glect some static or dynamic correlation. 

2. The symmetry dilemma. This refers to the fact 
that typical {e.g., Kohn-Sham (KS)) DFAs tend 
to break spin symmetry in strongly correlated sys¬ 
tems [H] whereas MR methods do not. One has 
therefore to choose between having physically pos¬ 
sible spin densities but unphysical energies, or im¬ 
proved energies but unphysical spin densities. Al¬ 
though one may opt on retaining the symmetries 
of the Hamiltonian, this choice results in massive 
static correlation errors for common DFAs [15] . 
Approaches to circumvent this issue include the use 
of functionals that take as inputs not only the den¬ 
sity but also the local pair density Hg ESHMl iQi 
m], or employing alternative spin densities defined 
by a transformation of the occupation numbers of 
the charge density matrix (M] |3ll [SS] |38l |42| . The 
former approach is formally justified by the work 
of Perdew et al. [IS], but can be computationally 
disadvantageous because computing the pair den¬ 
sity usually requires knowledge of the two-particle 
density matrix, which can be expensive to evalu¬ 
ate ISIIMIEI]- The latter approach is an ad hoc 
solution to the problem of computing the pair den¬ 
sity, but suffers from ambiguity |42l |43| in the pos¬ 
sible definitions of the alternative densities and is 
of course less rigorous. 

3. Problems of the MR method. Typical MR tech¬ 
niques of quantum chemistry have serious limita¬ 
tions. High computational cost and the lack of 
size-consistency and size-extensivity are probably 
the most prominent of these. Of the three main 
problems of MR-I-DFT, this is the most difficult to 
solve as it requires the development of novel, more 
efficient methods for handling static correlation. 
Alternative MR methods that have recently been 
proposed for use in MR-I-DFT include constrained- 
pairing mean-field theory |38l |50] (CPMFT), pro¬ 
jected Hartree-Fock [HI |321 ISI] (PHF), pair cou¬ 
pled cluster doubles [SHTTl [551 1^ (pCCD), and 
the density matrix renormalization group [231152] 
(DMRG). While they have afforded some encour¬ 
aging results, none of these alternatives is flawless: 
CPMFT-I-DFT often reduces (energetically) to un¬ 
restricted (U)HF-I-DFT unless DFT exchange is in¬ 
cluded in the mixture |38| (which is undesirable as 
it introduces self-interaction); PHF lacks size con¬ 
sistency and extensivity M; pCCD may not al¬ 


ways provide a complete description of the static 
correlation {e.g., when dissociating N 2 it goes to 
higher energy limit than UHF or CPMFT |1IJ[53]); 
and DMRG still bears some of the problems of 
traditional MR techniques {e.g., requires defini¬ 
tion of an active space, and size-consistency can 
only be guaranteed by enlarging this space). We 
should also mention here that attempts have also 
been made to add dynamic correlation to these 
techniques using approximations other than DFAs 
(see, e.g.. Refs. mUilHSI]), but, save a few excep¬ 
tions |5H1 [55] , the increase in cost is non-negligible 
in these approaches (and the perturbative ones can 
also introduce instabilities 0160]). 

Here, we propose a different way to elude these prob¬ 
lems of MR-I-DFT via CCO-fDFT combinations. The 
double counting is avoided by using the fact that CCDO 
and EDO have only correlations involving pairs of elec¬ 
trons with opposite spin; we show that the missing con¬ 
tributions can be added in terms of a multiple of the 
parallel-spin correlation from a density functional. The 
fact that we compute solely equal-spin correlation with 
the DFA allows us then to avoid the typical problems 
caused by the symmetry dilemma without resorting to 
symmetry breaking. Furthermore, by using an appropri¬ 
ate meta-GGA (generalized gradient approximation) for 
the DFA, the exactness of EDO for two-electron systems 
is maintained in BDO-f-DFT. Lastly, GGDO and EDO are 
size-consistent and size-extensive, and have polynomial 
cost, rather than the combinatorial cost of typical multi¬ 
configuration techniques. Extensions incorporating RPA 
correlation are also derived and studied. 


II. THEORY AND METHODS 

A. Singlet-Paired Coupled Cluster: CCDO and 
EDO 

Like other GG methods, singlet-paired coupled cluster 
doubles (GGDO) starts from an exponential ansatz |5] 

^[ 0 ] 

I'l'ccDo) = e 2 |$rhf), (2) 

where |$rhf) is a restricted Hartree-Fock reference de¬ 

terminant and the operator contains only the singlet- 
paired component of the double excitations; i.e., if we let 
i and a be indices for occupied and virtual orbitals, re¬ 
spectively, is 

= ( 3 ) 

ijab 

with 

Pij = (^) 


3 


To shed light on the idea behind CCDO (and this 
is helpful for understanding the justification for our 
CCDO+DFT combinations), it is instructive to consider 
standard CCD. The full double excitations cluster oper¬ 
ator T 2 used in the latter can be expressed as 


Tz +tW, 

(5) 

where is dehned above and 


ijab 

(6) 

where Qq may be written as a vector 
whose components are 

Qij = 

Qtj ~ Qij ~ 

(7) 

Q% — (^7^4 — 

(8) 

= ^(ffiTC4 + C4C4). 



results from CCSDO are not signihcantly different from 
those of EDO. 

The CCDO and EDO methods have certain advantages 
over traditional MR techniques. Eoth algorithms scale 
as the cost of CCDO is identical to that of CCD 

and the same holds for EDO and ED. While not precisely 
inexpensive, it is certainly better than the combinatorial 
scaling of typical MR techniques and, as CC-based meth¬ 
ods, CCDO and EDO may be able to take advantage of 
recently developed approaches [70] to decrease the cost 
of coupled cluster calculations. In addition, the use of 
exponential wavefunctions (Eqs. [^and|^ guarantees size 
extensivity. Size consistency is also satisfied by the ex¬ 
ponential ansatz provided that the reference determinant 
be size consistent. Lastly, it is easy to modify existing 
CCD and ED subroutines to implement CCDO and EDO, 
respectively. Easically, most restricted CCD codes use 
as fundamental variables the amplitudes of the full 
double excitations cluster operator T 2 ^ which may be ex¬ 
pressed as [63] 

^2 = ^ H cLci/jCj/SCic, (10) 

ijab q ;/3 


Note that acting on an empty state |) gives a singlet 
wavefunction (hence the name of singlet-paired coupled 
cluster), whereas Qjj |) yields pure triplet components 
with different m. We see then that CCD incorporates 
contributions from both the singlet- and triplet-paired 
components of T 2 . Thus, more correlation is taken into 
account in CCD as compared to CCDO. However, it is 
the combination of the singlet- and triplet-paired com¬ 
ponents what causes the failure of CCD (and other com¬ 
mon CC approximations) in strongly correlated systems, 
which results in the correlation becoming too large or 
even complex m- Ey removing from T 2 , CCDO re¬ 
linquishes a fraction of the correlation in exchange for 
safeguard against this breakdown. One may also choose 
to retain only the triplet-paired component in order to 
avoid the breakdown, but this alternative recovers less 
correlation than CCDO [5]. 

Just as CCD can be improved by using (approximate) 
Erueckner orbitals in the Erueckner doubles method [65l - 
[57] (ED), CCDO has an analogous extension in singlet- 
paired ED (EDO). Recall that the exact Erueckner 
orbitals are those that make the coefficients of singly- 
substituted determinants in the full configuration inter¬ 
action (FCI) wavefunction vanish [55]. The EDO ansatz 
is thus very similar to that of Eq. 

|^BDo)=e^^ 1‘1 ’bd), (9) 

except that the RHF reference is replaced by a determi¬ 
nant with approximate Erueckner orbitals |$bd) which 
zero out the amplitudes of single substitutions in a sub¬ 
space of single and double substitutions. We remark 
that a singlet-paired CCSD method (CCSDO) is possible 
too m, but (just as it happens for ED and CCSD [55] 1 


and, because the singlet-paired (triplet-paired) compo¬ 
nent of T 2 is symmetric (antisymmetric) with respect to 
the interchange of i and j, the afj' amplitudes (Eq. ^ of 
the singlet-paired component obey the following relations 


^ab ^ba ^ab ^ba 
^ij — — ^ji — ^ji 
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( 11 ) 


which implies that we can implement CCDO by replacing 
^i^jt iteration in a CCD 

code. Analogous remarks apply for the implementation 
of EDO starting from ED subroutines. 

The above discussion regarding CCDO and EDO suf¬ 
fices for the purposes of this paper. The interested reader 
may found further details about singlet-paired coupled 
cluster in Ref. |2]. 


B. Adding DFT Correlation to CCO Without 
Double Counting: CCO+DFT 

To combine CCDO and EDO with density functionals 
while avoiding double counting, we hrst note that the 
CCDO wavefunction can be expanded as 

roi 

1‘1’rhf) = 1‘1’rhf) + 72° 1‘1’rhf) H— ^2!^ I*^rhf) + • • ■ 
= |$rhf) H-, 

ijab 

( 12 ) 


and that the CCDO correlation energy is 


^CCDO 








ijab 


(13) 
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where = {ij\ab) is a two-electron integral in the Dirac 
notation. It is clear from these expressions that CCDO 
neglects contributions to the correlation from equal-spin 
excitations: only integrals involving pairs of opposite- 
spin electrons in occupied or virtual orbitals are weighted 
in the correlation energy. The same applies to EDO ex¬ 
cept, of course, that the orbitals i, j, a, and b would 
be (approximate) Brueckner, rather than RHF. One can 
thus think of adding equal spin correlation to CCDO using 
a density functional; the correlation of the parallel-spin 
density given by a functional contains, in principle, con¬ 
tributions from all excitations involving same-spin elec¬ 
trons only. Hence, we may, in the spirit of Eq. write 

^CCDO+pDFT ^ ^CCDO + ^ 

(14) 

where the “p” in pDFT is for “parallel-spin” and 
is the correlation of the spin-a density computed by a 
DFA. To evaluate Ef^^, we can use the exchange-like 
spin resolution of the correlation by Stoll et al. m 

(15) 

Also, because CCDO and EDO do not break spin symme¬ 
try and we will be working only with closed shells, the 
total CCDO-l-pDFT energy simplifies to 

2E^^^[n^, 0], (16) 

and likewise for EDO. 

Now, looking back at the triplet-paired component of 
T 2 that is excluded in CCDO (Eqs. 00 an d§ , we realize 
that, apart from the same-spin correlation, CCDO misses 
opposite-spin contributions from the m = 0 component 
of in Eq. 1^. However, for a closed shell, the 

m = -|-1, m = 0, and m = —1 channels of contribute 
all equally to the energy. In other words, the spin-up 
correlation associated with the m = -\-l element of 
is, by symmetry, identical to the opposite-spin correlation 
of the TO = 0 part. We may therefore incorporate the 
opposite-spin energy that is missing in CCDO-fpDFT by 
adding once more to the total energy 

-h [n^, 0], (17) 

where the “t” in tDFT indicates that the full contribu¬ 
tions from the triplet-paired component of T 2 are being 
taken into account. Again, the above analysis applies to 
EDO too. 

In this work, we employ Eqs. [TB] and [17| in a non-self- 
consistent manner: after a self-consistent CCDO (EDO) 
calculation, the equal spin correlation is obtained in a 
single-shot evaluation of a DFA using the densities from 
the reference RHF (Brueckner) determinant. This is a 
reasonable assumption because, when adding DFA dy¬ 
namic correlation to MR wavefunctions, the error in 
the approximate functional is often larger than the er¬ 
ror due to lack of self-consistency [72] • We benchmark 
both CCO-l-pDFT and CCO-l-tDFT with CCO referring 


to both CCDO and EDO. The motivation for testing 
both CCO+pDFT and CCO-ftDFT is that if the afj 
amplitudes are close to those of the full configuration- 
interaction (FCI), and the same-spin DFA correlation is 
fairly accurate, then we would anticipate CCO-ftDFT to 
outperform CC0-|-pDFT. However, if it is not the case 
that the afj amplitudes of CCO are similar to the exact 
ones, then CC0-|-tDFT may overcorrelate and not neces¬ 
sarily be better than CCO-fpDFT. 


C. Eluding the Symmetry Dilemma and Selecting 
an Adequate Density Functional 


So far we have identified the advantages of CCO 
over traditional MR techniques, and the derivation 
of CCO-I-DFT above explains how double counting is 
avoided; we are just left to deal with the symmetry 
dilemma. To illustrate how this issue can be avoided, 
consider the paradigm for static correlation in quantum 
chemistry: the H 2 singlet molecule at the dissociation 
limit |49| . For the exact, symmetry-adapted, wavefunc- 
tion, and are distributed equally in each atom, but 
the probability of finding both electrons at the same atom 
is zero. That is, the correlation is purely static and there 
is no dynamic correlation because the electrons are never 
close to each other. However, with symmetry-adapted 
densities, common DFAs yield a nonzero E^^^ that is 
caused (largely or completely) by the fact that semilocal 
approximations depending on the density and its gradi¬ 
ents “see” the densities of two different electrons on top 
of each other; they do not know that the local pair den¬ 
sity is zero everywhere. This problem is normally allevi¬ 
ated by breaking spin symmetry in a way that localizes 
an electron of a certain spin around a specific nucleus— 
so that the functional won’t see the densities of different 
electrons overlapping—or by introducing functionals that 
depend not only on the density, but also on the local pair 
density [48] . 

Because we evaluate the residual correlation in terms 


of parallel-spin correlation only (see Eqs. 16 and 171, 


the spurious opposite-spin correlation due to the use 
of symmetry-adapted densities when describing the dis¬ 
sociation limit of H 2 (and stretched bonds in general) 
disappears in CCO-fDFT. Therefore, there is no need 
for symmetry breaking or using the local pair density. 
Nonetheless, we should still consider the problem of self¬ 
interaction in approximate functionals, which produces 
artificial self-correlation and is present in the local den¬ 
sity approximation (LDA) and generalized gradient ap¬ 
proximations (GGAs). One can get rid of this issue, at 
least in the case of two-electron systems, by using meta- 
GGA functionals: DFAs that depend on the density, its 
gradient, and the kinetic energy density. We have there¬ 
fore chosen to test here combinations of GGO with two 
meta-GGAs: the Tao-Perdew-Staroverov-Scuseria m 
(TPSS) functional, and the recently developed strongly 
constrained and appropriately normed (SGAN) func- 
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tional of Sun et al. m- These functionals were selected 
based on the facts that (1) they are nonempirical; (2) 
they are free of one-electron self-interaction error, so that 
they yield no same-spin correlation for two-electron sin¬ 
glets and both BDO-fTPSS and BDO-I-SCAN will be ex¬ 
act for these systems; and (3) that despite these similar¬ 
ities, TPSS and SCAN are quite different because they 
were designed based on different paradigms. Specifically, 
TPSS respects the paradigms of describing correctly one- 
or two-electron densities and slowly varying densities m, 
whereas SCAN obeys all the 17 known exact constrains 
that a meta-GCA can obey, and is exact or near exact 
for a set of rare-gas atoms and certain nonbonded inter¬ 
actions [74] (the so-called “appropriate norms”). 


D. Improving the Spin Resolution of the SCAN 
correlation 

In Section 2.2 we used the spin resolution of by 

Stoll et al. [TT] , 0]) to extract 

the parallel-spin correlation. However, this decomposi¬ 
tion is correct only for fully spin-polarized densities and 
in the high-density limit of the uniform electron gas US). 
In other regimes of the uniform gas, it exaggerates the 
equal-spin correlation. It is possible, however, because of 
how the correlation is composed in SCAN, to formulate 
a better educated guess for the spin resolution of this 
functional. 

The dependence on the kinetic energy density, r(r) = 
^°“(1/2)|V ^/>i(r)p, is introduced in SCAN via the di¬ 
mensionless variable a m- 

a = (r-r'^)/T™'f>0, (18) 

where r'^ = |Vnp/8n is the single-orbital limit of r, 
and is the uniform density limit. SCAN constructs 
the correlation energy density £c by interpolation and 
extrapolation from the a = 0 (e°) and a = 1 (sj) limits. 
It is therefore natural to build its spin resolution in the 
same way: 


_ gin 


+ Ua) 


(19) 


the fraction of opposite-spin correlation density in the 
uniform density limit. Gori-Giorgi and Perdew m have 
worked out the spin resolution for the uniform electron 
gas: they determined fractions E„^i such that e'^'^ = 
ScFca' in terms of the local Wigner-Seitz radius Vg = 
(47rn/3)“^/^ and the spin polarization ^ = (n^ — n^)/n. 
Therefore, we calculate F^^{rsX) using Equation 9 of 
Ref. [75] , which is an interpolation between exact results 
for the high {vg —>■ 0) and low {vg —>■ oo) density limits 
that agrees with available quantum Monte Carlo data m 
in the range 0.8 < rs < 10 with C = 0. 

Note that integrating Eq. yields all the SCAN 
correlation that needs to be added in CC0-|-pDFT; for 
CCO-l-tDFT, this energy should be multiplied by 3/2 (see 
Section 2.2). Also, Eq. 20 gives no correlation for two- 
electron singlets because, for these systems, a = 0 and 
/c(0) = 1. Thus, BDO-fSCAN remains exact for two- 
electron systems when the parallel-spin correlation is re¬ 
solved using Eq. 


E. Modeling the Long-Range Correlation with the 
Random Phase Approximation 

Adding residual correlation to CCO via semilocal DFAs 
has the following pitfall: Semilocal functionals can’t cap¬ 
ture the long-range part of the correlation in van der 
Waals interactions [771 EH). Therefore, the description of 
these forces in CCO-l-DFT is essentially the same as that 
of CCO. Because van der Waals forces are very sensi¬ 
tive to dynamic correlation, and CCO misses a substan¬ 
tial part of it (the triplet pairing channel, to be precise), 
CCO-fDFT may be inadequate for simulating this sort of 
interactions. 

This problem can be solved by the technique of range- 
separation [791 - [M| : The interelectron Coulomb operator 
r^2 i® separated into a short-range (SR) component and 
its long-range (LR) complement 

1 I-erf(^ri2) ^ erf(^ri2) 

ri2 ri2 ri2 

SR LR 


where fda) (which is given in the Supporting Informa¬ 
tion of Ref. [71]) satisfies /c(0) = I, /c(I) = 0, and 
/c(oo) = —0.70, and likewise for the spin-down com¬ 
ponent. We note that = 0, since there is 

no parallel-spin correlation for two electrons in the same 
spatial orbital. Thus, the fraction of same-spin correla¬ 
tion in SCAN depends only on the fraction of same-spin 
correlation in the uniform density limit, -I- = 

(1 — F.fi)el. The total parallel-spin correlation energy 
density thus becomes 

ed+ed = il-fda)){l-Fu)el. ( 20 ) 

All of the terms in the above equation are already de¬ 
termined in the original SCAN functional except for 


where erf is the error function and /i a parameter defining 
the range separation. (Other partitions of are pos¬ 
sible, but the above choice is most convenient as it has 
been used to parametrize screened functionals and also 
facilitates the evaluation of two-electron integrals.) We 
can then evaluate the DFA correlation with the SR in¬ 
teraction only; the LR complement should be computed 
with an approximation capable of describing the long- 
range part of the van der Waals forces in a way that does 
not add double counting to CCO (and, desirably, that 
won’t increase the cost). As we shall explain shortly, a 
suitable approach that satisfies all of these requirements 
consists in applying the direct random phase approxima¬ 
tion (dRPA) in a manner similar to that by Janesko et 

al. [55] . 









6 


For the sake of clarity, let us consider how the RPA 
models correlation. The RPA requires the solution of an 
eigenvalue equation of the form [31 [S^ 153] 



where A, B, X, and Y are all of dimension ov x ov, 
with o and v being the number of occupied and virtual 
spinorbitals, respectively. In the full RPA, the elements 
of A and B are 

^ia.jb ~ i^a “t” (^3) 


Bia,jb = {ij\\ab). 


(24) 


Equations 22 -24 are used for the calculation of excited 
states. To obtain the ground state correlation energy, 
one notes that the Tamm-Dancoff approximation (also 
known as configuration interaction singles) also computes 
excited states with Eq. 22 but setting B = 0 (ie., it 


solves for AZ = Zu). While this approximation con¬ 
tains only excitation operators, the RPA also includes 
de-excitation operators that can be thought of as cor¬ 
relating the ground state. Thus, the RPA correlation 
energy for the ground state is normally written as mm 


E. 


RPA 


= -Tr(a; — A). 


(25) 


In the dRPA, the correlation is given by the same ex¬ 
pression, except the the exchange contributions are ne¬ 
glected (SUHHIHS] (he., {pq\\rs) is replaced by {pq\rs) in 
Eqs. 23 and|24|). The dRPA has been most useful for 


incorporating correlation in methods that already con¬ 
tain exchange (see, e.g., Refs. laisiiiiT] and references 
therein). Direct RPA also has the enormous advantage 
that the correlation is guaranteed to be real if the or¬ 
bitals obey the aufbau principle, which is not true for the 
full RPA [88]. As a matter of fact, the full RPA correla¬ 
tion becomes complex in the presence of an RHF insta¬ 
bility [Sni [3D| (negative eigenvalue in the Hessian). The 
dRPA does not suffer any of the instability problems dis¬ 
cussed above for full RPA or CCD because of the neglect 
of exchange. This minimalist description of the RPA cor¬ 
relation is sufficient for our purposes here; the interested 
reader is referred to reviews on the subject [86l [91] for 
further details. 

From this explanation, and the analysis in Section 2.2, 
it is straightforward to see how long-range dRPA correla¬ 
tion may be added to CCO without double counting. The 
correlation is still expressed by Eq. [^ only the elements 
of A and B need to be altered: 


■^iajb — (Ca ^i)bijSab 4 “ (^^1 ^ee 1^.)) ’ 


(26) 


Bia,jb = ( 27 ) 

where = 1 if the spin functions of the spinorbitals 
Xi(x) and Xj(x) are identical and = 0 otherwise. 


and (ijjugg|a6) indicates that the two-electron integral 
be evaluated with the long-range interaction 


(D>ee|a^) = / dxidx2X*(xi)x*(x2)^^(^^^^(^x 

J ri2 

Xa(xi)Xb(x 2 ). (28) 


By calculating the dRPA correlation with the A and B 
matrices as defined in Eqs. [M] and [27] only parallel spin 
correlation that does not overlap with the CCO correla¬ 
tion is obtained. This is the same strategy that we uti¬ 
lized in Section 2.2 to avoid the double counting between 
DFT and CCO. Therefore, the same considerations dis¬ 
cussed in Section 2.2 apply; one can add to CCO only the 
parallel-spin correlation, or the full triplet-paired com¬ 
ponent contributions by multiplying the total equal-spin 
correlation by 3/2. 

The RPA correlation describes correctly dispersion and 
van der Waals interactions [33] |13] , and is exact for long- 
range correlations [94]. In addition, Scuseria et al. [3] 
have shown that solving the RPA eigenvalue problem of 
Eq. 22 is equivalent to solving for T = YX“^ in a Riccati 
CCD equation 


B + AT + TA + TBT = 0, 


(29) 


and that in Eq. 25 can also be expressed as 


E. 


iRPA 


= oTr(BT). 


(30) 


Using the Cholesky decomposition of A and B, Equation 


29 


can be solved in 0{N'^) computational effort [3] | 
where N is the number of basis functions. This does not 
exceed the cost of CCDO. Hence, the dRPA correlation 
fulfills all of the requirements that we were looking for: 
it can be added to CCDO without double counting or 
increase in scaling, and it describes properly long-range 
interactions. As we do for E^^^, the dRPA correlation is 
obtained from a single-shot, post-SCF, calculation using 
the RHF or BD orbitals and added to the CCO energy. 

As mentioned above, E^^^ needs to be evaluated with 
the SR interaction only, in order to avoid double count¬ 
ing with the dRPA correlation. To the best of our knowl¬ 
edge, there are no parametrizations for the short-range 
correlation of TPSS or SCAN. However, codes for the 
range-separated parametrization of the LDA are avail¬ 
able m- Therefore, we use an approximate local scaling 
to estimate the SR meta-GGA same-spin correlation 


= 


i^{r) 




Vn,r)(i'^r, 

(31) 


where and are short- and full-range DFA 

correlation energy densities, respectively, taking spin 
densities as inputs. A similar LDA-based scaling has 
been shown to produce reasonable results in our recent 
work on range-separated hybrids of pCCD and density 
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TABLE I: Summary of the CCO+DFT methods tested here. The notation, closed-shell energy formulas, and relevant equations 
are given; CCO can refer to CCDO or BDO and the densities are from the RHF or Brueckner reference determinants, respectively; 
the “p” in pDFT is for parallel spin; the “t” in tDFT is for triplet-pairing component; is the spin-up short-range 

SCAN correlation using the spin resolution of Section 2.4; is the spin-up long-range dRPA correlation; and the rest 

of the notation is given in the text. 

Method Energy Formula Relevant Eqs. 

CCO+B?"''^ + See Sec. 2.6. 


CCO+pDFT 

^CCO 

+ 2B?“[nt,0] 

14 

15 


CCO-t-tDFT 

^cco 

+ 3B?“[nt,0] 

E 

F 


CCO+prSCAN 

^cco 

+ 2B'^^*”[nt,n4.] 

U 

20 


CCO-t-trSCAN 

^cco 


F 

20 


LC-CCO+pDFT 

^cco 

-P 0] + 

14 

F 

25-31 

LC-CCO+tDFT 

^cco 

-P 3By-°“[nt. 0] -P 


F 

25-31 

LC-CCO+prSCAN 

^cco 

+ rap] + 

14 

F 

25-31 

LC-CCO+trSCAN 

^cco 

+ 3By.f]:®‘^*”[rap, rap] + 

F 

F 

25-31 


functionals [5S]. Note that we have used the decompo¬ 
sition by Stoll et al. [TT] (Eq. 15) for determining the 
effect of considering only parallel spin correlation on the 
scaling factor; the Gori-Giorgi-Perdew m spin resolu¬ 
tion can’t be used for this because the fractions (see 
Section 2.4) are independent of /i. Nonetheless, 
may still be evaluated with the spin resolution of Eq. 
for the SCAN functional. Also, both 0) and 

are expected to exaggerate equal-spin corre¬ 
lation, resulting in error canc ella tion on the ratio that de¬ 
fines the scaling. Lastly, Eq. 


31 


maintains the = 0 
condition for two-electron singlets as long as be 

evaluated with an adequate meta-GGA. 


The separation of the electron-electron interaction 
with Eq.j^also requires one to define a range-separation 
parameter /i. In the case of standard LC-KS-DET, where 
Eq.j^is used to divide the exchange interaction only, the 
optimal /i is highly system and property dependent isg- 
uno]. However, here we are interested in splitting the 
correlation in SR and LR, and in this case fi is more 
universal: Eromager et al. |ini[ 1102] have demonstrated, 
based on physical arguments and numerical experiments, 
that the optimal fi for evaluating the SR correlation with 
semilocal DEAs is in the vicinity of /i « 0.4 au. There¬ 
fore, we set /r = 0.4 au in all of our calculations (note 
that this value is widely used in LC-KS-DET function¬ 
als such as, e.g., LC-wPBE [53]). Additionally, by fix¬ 
ing size consistency and extensivity are preserved; 
system-dependent definitions of gi not all can guarantee 
this [innKimj. We note that it is also possible to de¬ 
fine a physically-motivated local /i(r), but this has the 
caveat that locally-range separated hybrids are difficult 
and expensive to evaluate |104l 1105] . 

Lastly, we remark that range separation using long- 
range dRPA correlation is one way to introduce the 
long-range part of the van der Waals interaction in 
CCO-I-DFT; other approaches are possible. For example, 
specialized van der Waals functionals of varying degree of 
empiricism have been developed by various groups (see. 


e.g., Refs. |106H109] and references therein). These can 
be used in CCO-I-DFT as long as meaningful spin resolu¬ 
tion for the correlation of the functional exist. 


F. Summary of methods and notation 


The possible CCO-I-DFT combinations presented here 
and their corresponding closed-shell energy expressions 
are summarized in Table|T](in this work, we deal only with 
closed shells). All these methods have been implemented 
in a development version of Gaussian m . A vari¬ 
ant called CCO-I-E')?^"'" is also introduced where the full 
is added to the CCO energy. The purpose 
of including CCO-I-A)?^"’’ is to assess the effects of double 
counting, and how well do improved CCO-I-DFT combi¬ 
nations eliminate this problem. The explanation of the 
notation for the rest of the methods is as follows: A “p” 
before the functional name indicates that only parallel 
spin correlation from the DFA—and dRPA, for the meth¬ 
ods that include it—is added to CCO; the “t” that the full 
triplet-component of T 2 is added (see Section 2.2). The 
variant of SCAN called “rSCAN” uses the spin resolution 
from Eq.[^ rather than that of Eq.[^ Finally, the “LC- 
” prefix (which stands for “long-range corrected”) speci¬ 
fies that the short- and long-range correlation terms are 
evaluated with the DFA and the dRPA, respectively. The 
long-range correction used here should not be confused 
with the one used in standard LC-KS-DET; the former 
affects only the correlation, while the latter affects only 
the exchange. The exchange does not need corrections 
in CCO-I-DFT as it is calculated with the wavefunction 
method. We also remind the reader that the calculations 
are carried out in a non-self-consistent, post-CCO man¬ 
ner; a self-consistent implementation is possible but not 
the focus of this exploratory paper. 
























III. RESULTS AND DISCUSSION 
A. Description of short-range dynamic correlation 

For testing the CCO+DFT description of short-range 
dynamic correlation, we consider the equilibrium dis¬ 
tances (Re) and harmonic vibrational frequencies (uje) of 
a set of ten first- and second-row diatomics that has been 
studied in previous MR-I-DFT works [HI Hn]. These 
molecules are listed in Table |lT] they comprise well- 
known examples of single- and multiple-bonds, homonu- 
clear and heteronuclear diatomics; accurate experimental 
data are available from established databases inaiiia]. 
The motivation for studying Rg and We is as follows: 

increases in magnitude as the interatomic dis¬ 
tance is reduced, increasing the bond strength. Because 
is deliberately constructed to have double 
counting from this approach should predict too 

short bond lengths and too high frequencies. If the rest 
of the CCO-l-DFT combinations in Table |l] really work, 
then these problems would disappear and improvement 
over CCO should be observed. 

The bond lengths predicted by BDO and BDO-I-DFT 
are compared with accurate data in Table [TTl results by 
CCDO and CCDO-I-DFT are highly similar to those of 
their corresponding BDO counterparts. Likewise, LC- 
CCO-I-DFT and CCO-l-DFT are not significantly different 
in this case. Note that the H 2 molecule, for which BDO is 
exact, is included among the benchmark set. The accu¬ 
rate data are from experiments [11211113] and the differ¬ 
ence with BDO for H 2 (0.001 A) may be considered as an 
estimate of the effects of basis set incompleteness (and, in 
principle, also the Born-Oppenheimer approximation). 
As expected from the discussion above, BDO-l-if)?^"'" con¬ 
sistently underestimates the bond lengths (mean error, 
ME = —0.015 A), while the BDO-I-DFT methods that 
avoid double counting do not do this and furthermore 
improve upon BDO: the mean absolute errors (MAEs) 
of the different BDO-I-DFT combinations are 2-4 times 
smaller than the MAE of BDO (0.013 A), all providing 
very good results. Thus, the CCO-l-DFT methods ap¬ 
pear to be working as intended in the prediction of equi¬ 
librium distances. Perhaps the only unexpected result 
is that there does not appear to be much difference be¬ 
tween the different spin resolutions of SCAN (pSCAN 
and prSCAN), or between adding only parallel spin cor¬ 
relation or the full triplet component of T 2 (prSCAN and 
trSCAN). We also note that most of the largest errors 
occur in Li 2 . This diatomic is somewhat challenging to 
describe because, even at equilibrium, there may be some 
static correlation present |2T]. Nevertheless, the errors 
for the Re of Li 2 obtained here are much smaller than 
those reported for MCSCF (error = 0.258 A) and MC- 
SCF-^DFT (error « 0.150 A) in Ref. [HI- 

General observations are mostly similar for the har¬ 
monic vibrational frequencies shown in Table m Again, 
we exclude LC-CCO-f-DFT data because the effect of the 
long-range correction with dRPA correlation is negligi¬ 


ble on the calculated frequencies. A difference though 
is that there is more variation between BDO-I-DFT and 
CCDO-I-DFT, and hence we include data for both in this 
Table. The effect of double counting is as expected and 
CCO-l-if):’^'^ consistently overestimates the frequencies 
(ME = MAE « 70-80 cm'^ with TPSS, 40-60 cm"! with 
SCAN). Results from CCDO and BDO are close to each 
other (MAE = 22 cm“^ for both), but adding pTPSS 
to the former worsens results (MAE = 42 cm“^) while 
it does not for the latter (MAE = 25 cm“^). Adding 
pSCAN, prSCAN, or trSCAN to CCDO gives similar or 
slightly worse results than CCDO; for BDO, predictions 
become similar or better. The densities from the ap¬ 
proximate Brueckner orbitals may therefore be more ad¬ 
equate inputs for than the RHF densities. This 

observation may be related to the assertions about the 
physical relevance of Brueckner orbitals—connections to 
DFT |114| and even the Kohn-Sham orbitals |115| have 
been suggested in the literature. Note also that the ex¬ 
act Brueckner determinant can be defined as the Slater 
determinant having the largest overlap with FCI, and 
thus the densities can be presumed to be of better qual¬ 
ity than those from RHF. Overall, all of the combina¬ 
tions of CCDO or BDO with SCAN give satisfactory re¬ 
sults for both bond lengths and frequencies; mixtures 
with TPSS are also reliable for these properties except 
for CCDO-fpTPSS, which tends to overestimate the fre¬ 
quencies. 


As noted above, it is rather surprising that different 
spin resolutions, or whether one uses pDFT or tDFT, 
have little effect on the calculated bond lengths and 
frequencies. However, these properties depend on rel¬ 
ative energies only; total energies are certain to be af¬ 
fected by factors such as the incorporation of parallel 
spin vs. triplet pairing channel correlation (i.e., pDFT 
vs. tDFT). It is then logical to ask which of the differ¬ 
ent flavors of C CO- l-DFT provides more accurate total 
energies. Table IV compares the accurate [116] energy 
for the Neon atom with BDO and BDO-fDFT estimates; 
results from CCDO and CCDO-I-DFT differ by less than 
1 milihartree from their BDO analogues. The correlation 
missed by BDO in Ne is substantial (120 milihartrees). If 
raw E^^"^ [ri'f-, ni] is added, results are worsened due to 
overcorrelation (error « —225 milihartrees). In contrast, 
adding only parallel spin correlation improves the total 
energy, although underestimating it (errors = 41 and 
20 milihartrees for pTPSS and pSCAN, respectively). 
The best agreement with the accurate energy is obtained 
when the correlation from the full triplet pairing channel 
is added to BDO and using the spin resolution of Eq. [^ 
(errors = 2 and —10 milihartrees for tTPSS, and trSCAN, 
respectively), but there is overcorrelation if the spin res¬ 
olution of Stoll et al. m is used with SCAN (error = 
—29 milihartrees for tSCAN). These observations are in 
line with the theoretical arguments from Sections 2.2 and 
2.4; ie., that BDO misses the full triplet component of 
T 2 , and that the resolution by Stoll et al. m exagger¬ 
ates same-spin correlation in most regimes of the uniform 
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TABLE II: Accurate [11211113] equilibrium distances (_Re, in A) for a set of diatomics and the deviations (calculated — accurate, 
in miliangstroms) from these distances for BDO and BDO+DFT methods using a Cartesian 6-311++G(2df,2p) basis. ME is 
the mean error and MAE the mean absolute error. Results from CCDO and CCDO+DFT are very similar. 


Molecule 

Accurate 

BDO 




BD0+ 



^TPSS 

^SCAN 

pTPSS 

pSCAN 

prSCAN 

trSCAN 

H 2 

0.741 

1 

-2 

0 

1 

1 

1 

1 

LiH 

1.596 

6 

-10 

-12 

3 

-2 

1 

-5 

HF 

0.917 

1 

-10 

-5 

-5 

-1 

-1 

-2 

HCl 

1.275 

7 

-7 

-2 

0 

3 

4 

2 

Li2 

2.673 

-28 

-28 

-50 

-2 

-20 

-11 

-28 

C 2 

1.243 

4 

-9 

-6 

-1 

0 

0 

-2 

N 2 

1.098 

1 

-8 

-5 

-2 

-1 

-1 

-2 

F 2 

1.412 

15 

-34 

-26 

-14 

-6 

-2 

-11 

CI 2 

1.988 

41 

-19 

-9 

2 

12 

20 

9 

GIF 

1.628 

25 

-19 

-8 

0 

8 

12 

5 

ME 


7 

-15 

-12 

-2 

-1 

2 

-3 

MAE 


13 

15 

12 

3 

5 

5 

7 


TABLE III: Accurate [1121 I113[ harmonic vibrational frequencies (tUe, in cm“^) for a set of diatomics and the deviations 
(calculated — exact) from these values for CCO and CCO-I-DFT methods using a Cartesian 6-311-1—|-G(2df,2p) basis. ME is the 
mean error and MAE the mean absolute error. 

Molecule 

Accurate 

CCDO/BDO 



CCD0+/BD0+ 



^TPSS 

^SCAN 

pTPSS 

pSCAN 

prSCAN 

trSCAN 

H 2 

4401 

23/11 

65/60 

42/23 

23/11 

23/11 

23/11 

23/11 

LiH 

1406 

-6/-3 

32/30 

39/31 

0/4 

0/3 

-3/-11 

9/0 

HF 

4138 

36/19 

160/156 

103/84 

102/84 

22/24 

44/24 

58/38 

HCl 

2991 

-38/-38 

70/65 

35/21 

7/21 

-22/-21 

-32/-47 

-17/-29 

Li2 

351 

-3/32 

19/18 

11/10 

9/7 

9/26 

2/16 

12/5 

C 2 

1855 

22/12 

106/87 

77/58 

89/51 

40/29 

40/26 

50/36 

N 2 

2359 

36/19 

108/87 

77/54 

63/40 

37/25 

44/21 

51/29 

F 2 

917 

17/-31 

132/87 

110/65 

81/8 

64/9 

52/0 

76/22 

CI 2 

560 

-19/-19 

45/41 

32/27 

22/17 

7/10 

-2/-6 

10/5 

GIF 

786 

-24/-34 

54/41 

36/18 

20/6 

7/0 

-4/-15 

12/-4 

ME 


4/-3 

79/67 

56/39 

42/25 

19/11 

16/2 

28/11 

MAE 


22/22 

79/67 

56/39 

42/25 

23/16 

25/18 

32/18 


electron gas. 

Another test for dynamic correlation based on rela¬ 
tive energies, but more sensitive than the bond lengths 
and frequencies above, is the description of proton affini¬ 
ties. When a proton is added to a neutral molecule, the 
density is redistributed and dynamic correlation changes 
subtly. This effect is reasonably well described by den¬ 
sity functionals, which have average errors of about 1-3 
kcal/mol in small molecules [1171 1118] . Table com¬ 
pares the experimental (ZPE-corrected) proton affinities 
of eight molecules with those calculated by TPSS, SCAN, 
BDO, and BDO-t-DFT. The geometries (which are opti¬ 
mized at the MP2/6-31G(2df,p) level) and reference data 
were taken from Ref. |118| : the effects of geometric re¬ 
laxation are known to be negligible m- Again, we 
focus on BDO-I-DFT methods because CCDO-I-DFT re¬ 
sults are very similar and the long-range correction of 
LC-BDO-fDFT has little effect on the proton affinities. 
Pure BDO tends to overestimate the proton affinities (ME 


= 2.58 kcal/mol) and has the largest MAE of all methods 
(2.69 kcal/mol). TPSS results are similar (ME = MAE = 
2.65 kcal/mol), but SCAN is considerably better (MAE 
= 1.01 kcal/mol). Previous works noted that TPSS did 
not improve upon its Perdew-Kurth-Zupan-Blaha |119j 
predecessor in the prediction of proton affinities m- 
The fact that SCAN is better than TPSS at describ¬ 
ing the effect of protonization is reflected in BDO-I-DFT: 
the BDO-I-TPSS errors (MAE « 2.3 kcal/mol) are 2-3 
times larger than the BDO-I-SCAN errors. The best re¬ 
sults are obtained when the full triplet-pairing compo¬ 
nent of T 2 is added to BDO using the spin resolution of 
Section 2.4 (BD0-|-trSCAN, MAE = 0.66 kcal/mol). The 
observations here point toward a relationship between 
the accuracies of the DFA and BDO-I-DFA for a given 
property. Note, however, that not any functional can 
be utilized in CCO-t-DFT. For example, the Lee-Yang- 
Parr |120j (LYP) functional models all correlation as be¬ 
ing opposite-spin so that there is no parallel-spin 
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TABLE IV: Accurate |116| total energy for the Neon atom (Hartrees) and deviations (calcnlated — exact, in milihartrees) from 
this value for BDO and BDO+DFT methods using a Cartesian cc-pwCV5Z basis. Resnlts from CCDO and CCDO+DFT differ 
by less than 1 milihartree from their BDO counterparts. 







BDO-t 




Accurate 

BDO 

£;TPSS 

£.SCAN 

pTPSS 

tTPSS pSCAN 

tSCAN 

prSCAN 

trSCAN 

-128.938 

120 

-229 

-224 

41 

2 20 

-29 

33 

-10 


TABLE V: Accurate proton affinities (PA = E{M) — E{MH'^), in kcal/mol) for eight molecules and deviations (calculated — 
exact, in kcal/mol) from these valnes for TPSS, SCAN, BDO, and BDO+DFT using a Cartesian MG3S basis. The TPSS and 
accurate data are from Ref. m- Results from CCDO and CCDO+DFT are very similar to those from their BDO counterparts. 


BD0+ 


Molecule 

Accurate 

TPSS 

SCAN 

BDO 

pTPSS 

tTPSS 

pSCAN 

tSCAN 

prSCAN 

trSCAN 

NHa 

211.90 

1.70 

0.72 

2.84 

3.53 

3.88 

1.22 

0.41 

1.70 

1.13 

H 2 O 

171.80 

0.40 

0.24 

1.44 

2.43 

2.93 

0.39 

-0.14 

0.69 

0.31 

C 2 H 2 

156.60 

4.60 

3.83 

2.89 

4.01 

4.56 

0.95 

-0.02 

1.30 

0.50 

SiH4 

156.50 

3.10 

-0.34 

2.75 

0.46 

-0.68 

0.24 

-1.01 

0.84 

-0.12 

PHa 

193.10 

2.60 

-0.02 

4.55 

1.04 

-0.71 

-0.26 

-2.66 

1.09 

-0.65 

H 2 S 

173.70 

3.30 

1.45 

3.72 

2.74 

2.24 

1.04 

-0.30 

1.79 

0.83 

HCl 

137.10 

3.10 

1.50 

2.92 

3.18 

3.31 

1.42 

0.68 

1.83 

1.29 

H 2 

105.90 

2.40 

-0.01 

-0.44 

-0.44 

-0.44 

-0.44 

-0.44 

-0.44 

-0.44 

ME 


2.65 

0.92 

2.58 

2.12 

1.89 

0.57 

- 0.43 

1.11 

0.36 

MAE 


2.65 

1.01 

2.69 

2.23 

2.34 

0.75 

0.71 

1.22 

0.66 


correlation that can be extracted from LYP. 


B. Description of long-range dynamic correlation 

In all the benchmarks that we have hitherto discussed, 
the use of LC-CCO+DFT over CCO+DFT has negligible 
impact on the results. However, the effect of the long- 
range, dRPA-based, correction can be dramatic on the 
description of van der Waals interactions. One such ex¬ 
ample, the dissociation of a Helium dimer, is shown in 
Figure [2 The accurate data shown there has been taken 
from Ref. |122j . Standard Brueckner doubles (BD) cor¬ 
rectly describes the dissociation profile of He 2 . The corre¬ 
lation present in BD but not in BDO is crucial in this case; 
the latter is repulsive in the region near the correct mini¬ 
mum. The reason for the failure of BDO for long-range in¬ 
teractions can be explained as follows: BDO misses equal- 
spin correlation (plus an identical-in-magnitude m = 0 
component of which, in the short-range, is gener¬ 
ally much smaller than the opposite-spin contributions 
because same-spin electrons avoid each other due to the 
Pauli principle (antisymmetry of the wavefunction). In 
the long-range, however, electrons are far apart from each 
other and their spin is no longer relevant on their dy¬ 
namic correlation. Hence, the correlation missing in BDO 
becomes more important in the long range. 

From Figure[^we also see that BDO+TPSS adds a neg¬ 
ligible contribution to the BDO description of the He 2 dis¬ 
sociation; the same-spin TPSS correlation is practically 
zero along the whole curve. Surprisingly, BD+SCAN 



Internuclear distance (Bohr) 

FIG. 1: Counterpoise-corrected dissociation energy profiles 
for the Helium dimer calculated by various methods us¬ 
ing a Cartesian aug-cc-pVQZ basis. LC-BDO+trSCAN uses 
fj, = 0.4 au.; BDO+tSCAN has the same behavior as 
BDO+trSCAN. The curve marked as LC-BDO+trSCAN* uses 
the semiempirical prefactor of 1.5 on the dRPA correlation of 
Ref. |85| . The accurate data are from Ref. |122| . 


is radically different from BDO+TPSS, but fails badly 
in an opposite way: SCAN exaggerates the long-range 
parallel-spin correlation and leads to far too much bind¬ 
ing. This wrong behavior is corrected when the DFA 
correlation is evaluated in the short-range only and com¬ 
plemented with dRPA correlation in the long range. The 
LC-BDO+trSCAN and BD curves are close to each other. 
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TABLE VI: Accurate binding energy for the Neon dimer com¬ 
pared with counterpoise-corrected calculations by BDO and 
BDO-I-DFT methods at the experimental equilibrium distance 
of 5.841 an using a Cartesian aug-cc-pVQZ basis. The geome¬ 
tries and accurate binding energy were taken from Ref. |123 | . 
All energies are in kcal/mol. 


Method 

Binding Energy 

Accurate 

0.08 

BDO 

-0.03 

BDO+pTPSS 

-0.03 

BDO+tTPSS 

-0.02 

BDO+pSCAN 

0.03 

BDO+tSCAN 

0.06 

BDO+prSCAN 

0.01 

BDO+trSCAN 

0.03 

LC-BDO+pTPSS 

0.00 

LC-BDO+tTPSS 

0.02 

LC-BDO+prSCAN 

0.01 

LC-BDO+trSCAN 

0.03 

LC-BDO+tTPSS* 

0.05 

LC-BDO+trSCAN* 

0.05 


* Using the semiempirical prefactor of 1.5 on the 
dRPA correlation of determined on Ref. 185 II . 


indicating that adding the full triplet component of T 2 
via our SR-SCAN/LR-dRPA model correctly incorpo¬ 
rates the correlation that is missing in BDO. It is worth 
pointing out that Janesko et al. found that using 
an empirical scaling of 1.5 on the dRPA correlation im¬ 
proves the agreement of LC-LDA-bdRPA with accurate 
data for noble gas dimers, and that here too we observe 
improvement if this factor is used (LC-BDO-|-trSCAN* 
curve). Note also that we do not show data for all pos¬ 
sible method combinations in Figure but the perfor¬ 
mance of the techniques not shown can be inferred by 
the reader due to the fact that many BDO-I-DFT vari¬ 
ants differ only by a multiplicative factor on the correla¬ 
tion added to BDO. For example, the LC-BDO-|-prSCAN 
curve (not shown in Figure]^ would be qualitatively sim¬ 
ilar to the LC-BDO-btrSCAN one, but less accurate in 
quantitative terms. 

On Table |VI| the binding energies of the Neon dimer 
calculated by BD and BDO-I-DFT at the experimental 
equilibrium distance of 5.841 au are compared with ac¬ 
curate estimates from standard databases |123j . Once 
more, BDO fails to bind the noble gas dimer, and so 
do BDO-bpTPSS and BDO-btTPSS. Remarkably though, 
all combinations of BDO with SCAN bind Ne 2 , show¬ 
ing best agreement with the reference energy when the 
full triplet-pairing component of T 2 is included. When 
the long-range DFA correlation is replaced by dRPA cor¬ 
relation, variants using TPSS and SCAN all bind the 
dimer, but tSCAN-based methods are most accurate. 
LC-BDO-bpTPSS provides the smallest binding energy 
(0.004 kcal/mol) and LC-BD0-|-trSCAN the largest one 
(0.03 kcal/mol). As occurs for He 2 , if the semiempiri- 


TABLE VII: Accurate binding energies for ethylene (C 2 II 4 ) 
and acetylene (C 2 H 2 ) dimers compared with counterpoise- 
corrected calculations by BDO and BDO-I-DFT methods using 
a Cartesian MG3S basis. The geometries and accurate data 
were taken from Ref. |I23| . All energies are in kcal/mol. 


Method 

Binding Energy 

(C2H2)2 

(C2H4)2 

Accurate 

1.34 

1.42 

BDO 

0.38 

-0.42 

BDO+pTPSS 

0.64 

0.07 

BDO+tTPSS 

0.76 

0.31 

BDO+pSCAN 

1.13 

1.03 

BDO+tSCAN 

1.50 

1.76 

BDO+prSCAN 

0.95 

0.23 

BDO+trSCAN 

1.24 

0.56 


cal factor of 1.5 determined on Ref. [85] is used to scale 
^dRPA agreement between LC-BDO-btDFT and ex¬ 
periment is improved. 

It should also be mentioned that there are some types 
of weak interactions that are captured by semilocal func¬ 
tionals, as these can describe the intermediate-range part 
of the van der Waals interaction As an 

example, we consider the tt-tt interactions in ethylene 
and acetylene dimers with geometries from a standard 
dataset m- Accurate data for the binding energies of 
these species are compared with BDO and BDO-fDFT 
predictions on Table |VII| At the accurate equilibrium 
distance, BDO significantly underestimates the binding 
energies for both dimers and predicts only the acety¬ 
lene dimer to be bound. All of the BDO combinations 
with DFT improve the former’s results and bind both 
dimers, but the best agreement with the reference data is 
given by BDO-ftSCAN and BD0-|-trSCAN. It was noted 
on Ref. [74] that SCAN was more accurate than TPSS for 
describing weak interactions; a consequence, most likely, 
of the appropriate norming of SCAN. Hence, we observe 
again a correlation between the adequacy of the DFA for 
a given property, and that of BDO-fDFA methods for the 
same type of calculation. 

In this Section, we have focused on BDO-I-DFT re¬ 
sults; CCDO-l-DFT results have been omitted. The rea¬ 
son for this is that, in our experience, BDO-l-DFT and 
CCDO-l-DFT are largely similar when describing prob¬ 
lems dominated by dynamic correlation (with BDO being 
slightly better). Thus, the observations for CCDO-based 
methods are analogous to those of their BDO counterparts 
and discussion about the former has been left out to avoid 
repetitiveness and long-windedness. Nevertheless, as we 
discuss next, there can be important differences between 
CCDO and BDO when static correlation is involved. 















12 




FIG. 2: Dissociation energy profiles for boron hydride calculated by various methods using a Cartesian cc-pVTZ basis. 


C. Simultaneous description of static and dynamic 
correlation 

Molecular dissociations are one of the paradigms of 
static correlation in quantum chemistry. Figure shows 
the dissociation energy curves for boron hydride calcu¬ 
lated by various standard and CCO-brSCAN methods; 
CCO-I-TPSS and CCO-I-SCAN variants are similar but 
slightly higher and lower, respectively, in energy than 
CCO-l-rSCAN. The UCCSD(T) (unrestricted CCSD with 
perturbative triples) data in this Figure may be consid¬ 
ered as an accurate reference, as this technique is con¬ 
siderably parallel to FCI for this system |126) (nonparal- 
lelity error « 2 kcal/mol). Note, however, that restricted 
CCSD(T) experiences a breakdown at B-H distances 
greater than about 5.5 au. As noted in the previous sec¬ 
tions, CCDO- and BDO-based methods are highly similar 
when the problem is dominated by dynamic correlation 
(in this case, near equilibrium, although the existence of 
a UHF solution slightly lower in energy than RHF sug¬ 
gests the presence of some static correlation too; similar 
observations have been noted by Pollet et al. m)- How¬ 
ever, this is no longer true at large bond lengths, where 
CCDO becomes significantly higher in energy than BDO. 
In general, the BDO description of bond breaking is sig¬ 
nificantly better than that of CCDO. This is easily un¬ 
derstood by considering the dissociation of H 2 : BDO is 
exact but CCDO is not because, in order to be exact for 
two-electron systems, singles contributions must be in¬ 
cluded either explicitly or via the Brueckner orbitals. As 
the total correlation increases when the bond is extended 
beyond the Coulson-Fischer point, the singles contribu¬ 
tions that are missing in CCDO become more important. 
This analysis extends to other chemical bonds in general. 

Although the BDO-|-trSCAN curve in Fig. [^provides a 
very good description of the BH dissociation (and with¬ 
out resorting to symmetry breaking, like UCCSD(T)), 
LC-BDO-|-trSCAN exhibits wrong behavior at large bond 
lengths. The dRPA correlation increases too much be¬ 
yond i?B-H « 6 au, resulting in a “bump”, similar to 


that of RCCSD(T), characteristic of common perturba¬ 
tive and RPA-based approaches. In the case of LC- 
CCDO-btrSCAN, the substantial dRPA correlation ac¬ 
tually helps to correct, to a certain extent, the too- 
high dissociation limit of CCDO-l-trSCAN. When break¬ 
ing multiple bonds, coupled cluster methods for strong 
correlation—CCDO, BDO, and pCCD—all tend to a too- 
high energy limit (sometimes higher than UHF, e.g., in 
N 2 ; see Figure [^. Thus, adding dRPA correlation in a 
manner similar to the one done here may provide a route 
to alleviate this problem of singlet-paired coupled cluster 
techniques. This, however, requires further investigation 
that is beyond the scope of the present paper. In the case 
of the dissociation of N 2 shown in Figure there also 
appears to be an improvement of the too-high dissocia¬ 
tion energy limit of CCDO by the long-range dRPA cor¬ 
relation. LC-CCDO-btrSCAN is quite close energetically 
to UCCSD and binding energies are no longer exagger¬ 
ated. The BDO-fDFT curves are above LC-CCDO-I-DFT 
at large internuclear distances, but LC-BDO-I-DFT (not 
shown) leads to too much correlation at dissociation and 
a pronounced bump similar to that seen in boron hy¬ 
dride. Note also that although this is one of the cases in 
which BDO goes to a too high limit, CCSD and CCSDT 
fail very badly while BDO-I-DFT provides a reasonable 
description of the dissociation (and is also accurate near 
equilibrium, as evidenced by the data in Tables |TT] and 

g. 

Another problem for which traditional coupled cluster 
fails badly is the deformation of H 4 from a rectangu¬ 
lar to a square geometry [5]. Figured shows this model 
(which has been studied extensively [127H131] ! schemati¬ 
cally. The H atoms are confined to a circle of radius 3.284 
au and preserve D 2/1 symmetry, so that the geometry de¬ 
pends on a single parameter: the angle 0 that bisects 
two lines connecting opposing H atoms. This Figure also 
shows the dependence of the energy on 0 for various 
methods, including FCI data from Ref. |13I) . It is seen 
that CCSD predicts a trend that is opposite to that of 
FCI, with a minimum at 0 = 90° (square geometry). 
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Internuclear distance (Bohr) 



FIG. 3: Dissociation energy profiles for the N 2 molecule cal¬ 
culated by various methods using a Cartesian cc-pVDZ basis. 


where it overcorrelates the most. In contrast, EDO and 
BDO-I-DFT methods have good qualitative and quantita¬ 
tive agreement with FCI, albeit they exhibit a disconti¬ 
nuity in the first derivative of the energy with respect to 
0 at 0 = 90° that is not present in the FCI curve. Of the 
BDO-I-DFT methods, only BD0-|-rSCAN data are shown 
in Figure but using TPSS gives very similar results, 
while BDO-bSCAN yields slightly lower total energies. 

A benchmark for which KS-DFT methods fail catas¬ 
trophically is the Beryllium isoelectronic series (z.e., 
X(Z-4) + 

ions). In this system, the angular 
static correlation increases linearly with increasing nu¬ 
clear charge Z] this effect is poorly described by common 
DFAs (even if one tries to break spin symmetry, be cause 
the ions are RHF —>■ UHF stable). Table 


VIII 


compares the accurate energies |132j for this series with 
energies from TPSS, SCAN, BDO, and their combina¬ 
tions; CCDO-l-DFT and BDO-l-DFT results are not sig¬ 
nificantly different in this case. The nonparallelity errors 
(NPE)—the difference between the maximum and min¬ 
imum errors—are also provided as an estimate of how 
well is the linear trend of the correlation described. As 
expected, TPSS and SCAN fail badly to capture this 
trend with NPEs of 40.6 and 87.9 milihartrees, respec¬ 
tively. Much better results are given by BDO (NPE = 



r = 3.284 au 



FIG. 4: Geometry of H 4 on a circle with a radius of 3.284 au. 
Dependence of the total energy on the angle 0 for various 
methods using a Dunning DZP basis. The FGI data are from 
Ref. |I31| . 


5.2 milihartrees), although there is underestimation of 
the energy due to the missing same-spin correlation (ME 
= MAE = 7.1 milihartrees). All of the BDO-I-DFT com¬ 
binations capture correctly the dependence of the energy 
with respect to Z, and are similar or better than BDO 
with NPEs in the range of 2.0-7.3 milihartrees. The best 
results are given by BD0-|-pTPSS (ME, MAE, and NPE 
of —2.5, 2.5, and 1.2 milihartrees, respectively). The 
BDO-I-DFT methods tend to overestimate the correla¬ 
tion, which may be a result of the approximate func¬ 
tional and/or of the resolution of the equal-spin correla¬ 
tion. Indeed, the latter appears to be a factor because 
BDO-bSCAN overcorrelates more than BD0-|-rSCAN. It 
is also possible that the overcorrelation in -l-tDFT meth¬ 
ods for the ions be in part an artifact of the 

small size of the system; total energies for the isovalent 
Mg atom (Table [D^ are most accurate with BD0-|-tTPSS 
and BDO-ftrSCAN and trends are similar to those ob¬ 
served for the Ne atom on Table |IV[ with BDO-f DFT 
providing great improvement over BDO. 

In the tests for static correlation thus far analyzed, 
BDO is able to provide a qualitatively correct description 
of the problem at hand. A challenging benchmark for 
which this is not the case is the dissociation of the Beryl¬ 
lium dimer. This dissociation is chemically different from 
the breaking of typical bonds discussed above: Be 2 has 
a formal (MO-theory) bond order of zero; it is weakly 
bound (« 2 kcal/mol) but not by dispersive van der 
Waals forces, but by a mixture of dynamic and static cor- 
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TABLE VIII: Accurate | 132 | total energies (in Hartrees) for the Beryllium isoelectronic series and deviations (calculated — 
exact, in milihartrees) from these values for TPSS, SCAN, BDO, and BDO+DFT using a Cartesian cc-pCVQZ basis. Results 
from CCDO and CCDO+DFT are very similar to those from their BDO counterparts. 









BDO-I- 



Species 

Accurate 

TPSS 

SCAN 

BDO 

pTPSS 

tTPSS 

pSCAN 

tSCAN 

prSCAN trSCAN 

Be 

-14.6674 

-4.2 

17.3 

4.0 

-2.7 

-6.0 

6.0 

-11.1 

-3.7 

-7.5 

B+ 

-24.3498 

3.0 

31.8 

6.0 

-2.4 

-6.5 

-6.8 

-13.2 

-4.2 

-9.3 

C2+ 

-36.5352 

8.8 

45.2 

6.2 

-3.2 

- 7.9 

-8.4 

-15.7 

-5.8 

-11.7 

N=>+ 

-51.2234 

15.6 

59.9 

7.1 

-3.0 

-8.0 

-8.8 

-16.7 

-6.2 

-12.9 

o^+ 

-68.4128 

22.8 

75.4 

8.3 

-2.3 

- 7.5 

-8.5 

-17.0 

-6.1 

-13.3 

p5-t- 

-88.1022 

29.5 

90.2 

8.7 

-2.2 

- 7.7 

-8.9 

-17.7 

-6.6 

-14.2 

Ne'‘+ 

-110.2921 

36.4 

105.1 

9.2 

-2.0 

-7.6 

-9.0 

-18.1 

-6.8 

-14.8 

ME 


16.0 

60.7 

7.1 

-2.5 

- 7.3 

-8.1 

-15.6 

-5.6 

-12.0 

MAE 


17.2 

60.7 

7.1 

2.5 

7.3 

8.1 

15.6 

5.6 

12.0 

NPE 


40.6 

87.9 

5.2 

1.2 

2.0 

3.0 

7.1 

3.1 

7.3 

TABLE IX: Accurate [116j total energy for the Magnesium atom (Hartrees) and deviations (calculated 
from this value for BDO and BDO-I-DFT methods using a Cartesian cc-pwCV5Z basis. 

— exact, in 

milihartrees) 








BD0+ 




Accurate 

BDO 


pTPSS 

tTPSS 


pSCAN 

tSCAN 


prSCAN 

trSCAN 

-200.053 

143 


48 

1 


24 

-36 


38 

-15 


relation m- Because of this, common coupled cluster 
and perturbation methods fail for this system; multiref¬ 
erence wavefunctions may or may not work depending on 
the size of the active space and on how dynamic correla¬ 
tion is treated |133til35] . Figure]^ shows that BDO also 
fails badly for Be 2 ; the potential energy curve is repulsive 
near the accurate equilibrium distance (determined by an 
explicitly correlated MR wavefunction, which agrees with 
experiment |I34l 1135] ). In contrast, BDO-I-DFT methods 
bind the Be dimer, with BD0-|-tTPSS giving particularly 
good qualitative and quantitative agreement with the ac¬ 
curate reference data. In this case, BDO- and CCDO- 
based provide similar results. LC-BDO-|-tTPSS affords a 
good description of the dissociation when the semiempir- 
ical prefactor of 1.5 on the dRPA correlation of Ref. [85] 
is used; otherwise the relative energy at the equilibrium 
distance may be too high (albeit there is still substantial 
improvement over BDO). It has been noted that con¬ 
nected triples are needed to obtain sizable bonding |136j ; 
the ability of BDO-t-DFT to do this is thus a good in¬ 
dicative that the way in which the methods are mixed is 
adequate. The dissociation of Be 2 is an example of the 
importance of the correlation missing in CCO for obtain¬ 
ing qualitatively correct results in difficult cases where 
a simultaneous, balanced description of both static and 
dynamic correlation is necessary. 


IV. CONCLUSIONS, OUTLOOK, AND 
POSSIBLE IMPROVEMENTS 

We have presented techniques to incorporate the cor¬ 
relation that is missing in singlet-paired coupled cluster 
via density functionals or combinations of these with the 
dRPA. These methods are nonemprical and physically 
motivated, and the benchmarks studied here demon¬ 
strate that they are capable of describing static and 
dynamic correlation (including long-range weak interac¬ 
tions) without symmetry breaking. Typical problems of 
MR-l-DFT methods such as double counting and the sym¬ 
metry dilemma are avoided in CCO-I-DFT. In general, 
CCDO- and BDO-based methods provide similar results 
for weakly correlated systems, but BDO-t-DFT is prefer¬ 
able for problems with static correlation such as bond 
breaking. Similarly, both CCO-I-TPSS and CCO-I-SCAN 
yield good overall results, although SCAN is better at 
capturing long-range weak interactions. Addition of the 
full triplet-pairing component of T 2 to CCO (-l-tDFT) is 
frequently better than adding just the parallel-spin cor¬ 
relation (-|-pDFT), in agreement with theoretical argu¬ 
ments outlined in Section 2.2. 

It is often said of traditional coupled cluster meth¬ 
ods that they give “the right answer for the right rea¬ 
son.” Yet, these approximations fail in the presence of 
strong correlation. Here, LC-CCDO-fDFT, for example, 
can be considered as an attempt to avoid this failure 
while obtaining the right answer for the right reason via 
a different route: CCDO correctly describes most of the 
opposite spin correlation without breaking down in the 
strong correlation regime; DFAs are most accurate in the 
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Be 2 /aug-cc-pVTZ 



CCDO 

♦ 

LC-CCDO+trSCAN 

-•0-- 

CCDO+trSCAN 


LC-BDO+tTPSS* 

— 

CCDO+tTPSS 

— 

Accurate 



-i—I—I—I—I—|- 

5 6 7 

Internuclear distance (Bohr) 


FIG. 5: Dissociation energy profiles for the Beryllium dimer 
calculated by various methods using an uncontracted Carte¬ 
sian aug-cc-pVTZ basis. The accurate data are from explic¬ 
itly correlated ri 2 MRCI calculations from Ref. [135| . The 
CCSD data are from Ref. | 134| . The curve marked as LC- 
BD0-|-tTPSS* uses the semiempirical prefactor of 1.5 on the 
dRPA correlation of Ref. | 85| . 


short-range |137j , while RPA is most accurate in the long- 
range [92H94] ; and we add specihcally the correlation ab¬ 
sent in CCDO from these last two. Nonetheless, there 
is room for improvement over the techniques presented 
here: 

• The major drawback of CCDO and BDO is their 
0{N^) scaling, which determines the cost of the 
combinations presented here. Additionally, the 
DFA correlation added to CCDO/BDO is evalu¬ 
ated in terms of a multiple of the same-spin cor¬ 
relation. This has the disadvantage that resolu¬ 
tions for the parallel-spin correlation are imper¬ 
fect, and also that same-spin correlation typically 
makes up for a larger fraction of the total correla¬ 
tion in the uniform electron gas as compared to 
molecules [30l [71], which may result in too low 
total energies in CCO-I-DFT. Both of these prob¬ 
lems can be hxed by appropriate combinations of 
DFT with pCCD, rather than BDO and CCDO. 


The scaling of pCCD is only 0{N^) (neglecting 
the basis transformation of the two-body interac¬ 
tion my Although pCCD misses more dynamic 
correlation than CCDO, the intrapair correlation is 
described almost exactly; pCCD closely reproduces 
seniority-zero FCI—an optimal linear combination 
of all configurations that preserve electron pairs m 
{i.e., all seniority zero configurations; see Ref. |138j 
for a detailed explanation of the concept of senior¬ 
ity). Hence, following the philosophy of the present 
work, we can add the correlation missing in pCCD 
by switching the discussion from the singlet/triplet¬ 
pairing components of CCO to intra/inter-pair cor¬ 
relation. Ways of extracting intra/inter-pair cor¬ 
relation from DFAs have been discussed by Savin 
et al. HMj. Likewise, adding dRPA correlation to 
pCCD without double counting can be done with a 
similar strategy to the one done here by including 
only integrals that “break” electron pairs (seniority 
nonzero) in the B matrix of dRPA. 

• Here, the DFA correlation is added to CCO in a 
non-self-consistent manner. Although the effect of 
self-consistency is often small when adding dynamic 
DFA correlation to a wavefunction it is possi¬ 
ble that self-consistency can improve CCO-I-DFT. 
In particular, it would be possible to make the 
singlet- and triplet-pairing components of T 2 to 
“talk” to each other—and without introducing the 
risk of breakdown in strongly correlated systems— 
by including an effective one-body potential from 
the DFA in core Hamiltonian of CCO. 

• We have relied here on educated guesses for extract¬ 
ing the parallel spin and short-range correlation 
from existing functionals. Thus, our CCO-fDFT 
methods can be refined by parametrizing a parallel 
spin correlation functional complementary to CCO, 
and doing a rigorous parametrization of said func¬ 
tional for the range separation. Related to this, 
LC-CCO-I-DFT may also benefit from a hner tun¬ 
ing of the range separation parameter. 

Thus, it is likely that the already good results obtained 
here can be improved further, and that the cost of 
CCO-I-DFT be reduced by using pCCD-|-DFT combina¬ 
tions. Indeed, many other CCO-I-DFT (or CCO-fRPA) 
mixtures can be developed following the general strate¬ 
gies used here. Currently, we are working on some of 
these extensions and on testing CCO-I-DFT methods on 
larger systems; preliminary results have been encour¬ 
aging. We hope that the ideas and results presented 
here can stimulate the further development of CCO-I-DFT 
methods. 
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